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The XY model with quenched random disorder is studied by a zero temperature domain wall 
renormalization group method in 2D and 3-D. Instead of the usual phase representation we use the 
charge (vortex) representation to compute the domain wall, or defect, energy. For the gauge glass 
corresponding to the maximum disorder we reconfirm earlier predictions that there is no ordered 
phase in 2D but an ordered phase can exist in 3D at low temperature. However, our simulations 
yield spin stiffness exponents 6a ~ —0.36 in 2D and 9s « +0.31 in 31?, which are considerably 
larger than previous estimates and strongly suggest that the lower critical dimension is less than 
three. For the ±J XY spin glass in 3-D, we obtain a spin stiffness exponent 9s ~ +0.10 which 
supports the existence of spin glass order at finite temperature in contrast with previous estimates 
which obtain 9s < 0. Our method also allows us to study renormalization group fiows of both the 
coupling constant and the disorder strength with length scale L. Our results are consistent with 
recent analytic and numerical studies suggesting the absence of a re-entrant transition in 2D at low 
temperature. Some possible consequences and connections with real vortex systems are discussed. 
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I. INTRODUCTION 

The XY model with quenched random phase shifts 
as a model for a superconducting glass has been inten- 
sively investigated over the last decade, focusing on the 
so-called gauge glass model which corresponds to the case 
with maximal disorder. Since a transport current exerts 
a force on a flux lattice, it tends to move in response 
which causes dissipation of the current. The existence of 
disorder, which destroys the flux lattice structure, is es- 
sential to pin the vortices in order for a supeeconducting 
phase to exist in a high Tc superconductortlElH. Although 
there exists no proof whether or not the gauge glass and 
vortex glass are in the same universality class, it is of 
interest as the simplest model of a disordered supercon- 
ductor and is still not understood despite all the effort 

From numericafliSBifl and experimentalS studies, it 
is believed that the gauge glass has no ordered phase at 
any finite temperature in two dimensions. In three di- 
mensions, numericj]] domain wall renormalization group 
(DWRG) studieslllrEj indicate that the lower critical di- 
mension seems to be close to three. However the situa- 
tion is less conclusive, since the simulations are limited 
to small system sizes. Finite temperature Monte UsfJA 
studies yield a transition temperature Tc/ J ^ 0£LxMx3 
which is difficult to reconcile with DWRG studiesElffliJ as 
these studies imply that the lower critical dimension for 
superconducting glass order is close to three. Experimen- 
tally there is also some evidence for a finite tempecaiure 
phase transition to a superconducting glass phaseoEj. 

The Hamiltonian of the XY model with random 
quenched disorder can be written as 

H=J2v{e,~e,^A,,) (1) 



where V{(f>) is an even, 27r periodic function of (j) with 
a maximum at (j) — t: and minimum at (j) — 0, usually 
taken as V{(j)ij) — —JijCOs{(f>ij). The sum is over all near- 
est neighbor pairs of sites and the coupling constants, Jij , 
are assumed uniform, Jij = J > 0. The random bond 
variables Aij, which are responsible for the randomness 
and frustration, are taken to be independent and uni- 
formly distributed in {—air, air] with < a < 1. For a 
gauge glass, 9i is the phase of the superconducting order 
parameter at site i of a square lattice in 2D and a simple 
cubic lattice in 3-D. The random bond variables Aij are 
taken to correspond to maximal disorder with a = 1. An 
external field applied to an extreme type II superconduc- 
tor induces a uniform component A^j = (27r/<I>o) J!' A-dl 
where A is the vector potential of the applied field and 
<f>o = hc/2e is the quantum of flux. In this work, we 
take A'^j = 0, corresponding to zero applied field. Unless 
explicitly stated, we consider an unscreened system with 
a — 1 corresponding to maximal disorder. The Hamil- 
tonian of Eq. (|^) also describes the XY j-magnet with 
random Dzyaloshinski-Moriya interactional^ and alsOj-a 
Josephson ju]ae|tiQp-. array with positional disorder. EZlE£l 
These studiesESEZEa showed that the existence of weak 
disorder {a ^ 1) does not destroy an ordered phase 
at intermediate temperature but predict a re-entrant 
transition to a disordered phase at low tenTperatiye in 
two dimcH^ions. However, recent analyticE3E3cd'L3 and 
numcricalQ studies suggest the absence of a re-entrant 
transition and that there exists an ordered phase for 
T < Tc{a) when a < ac- 

When the random bond variables Aij are restricted 
to or TT with equal probabilities, this model reduces 
to the ±J XY spin glass, which is believed to be in 
a different universality class due to the additional re- 
flection symmetryo which is absent in the case of uni- 
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formly or Gaussian distributed Aij. An XY spin glass 
may have both spin and chiral glass order associated 
with rotational and reflection symmetry, respectively. 
It has been suggested that, in 2D and 31?, spin and 
chiral variablp-pxlemaiple at long distances and order 
independentlyE3Ej'L£lLZl, and the lower critical dimensions 
are di > A for spin glass order and < 3 for chiral glass 
order. However, the decoupling sceaario contradicts the 
analytic studies on a ladder latticeEa, on a tube latticed, 
and on a 2D lattice with a special choice of disordered. 
Recent numerical simulation£3 also suggests di for a spin 
glass order may be close to three. 

In this paper, we re-investigate the possibility of 
an ordered phase at small but finite temperature T 
by a mHHpHcal domain wall renormalization group 
(DWRG)lliHi3, or defect energy scaling. The domain wall 
or defect energy of the system is computed by using the 
Hamiltonian in the Coulomb gas (vortex) representation, 
which is more convenient for numerical work as it elimi- 
nates spin wave contributions to the energy. Although 
the conventional DWRG method can handle only the 
scaling of the coupling constant J{L) at scale L, which 
is proportional to the domain wall or defect energy, our 
method enables us to study the flows of both the cou- 
pling constant and the disorder strength, A{L), at length 
scale In. We apply this to the case of general disorder 
strength, < a < 1. The outline of the paper is as fol- 
lows. In Section || we discuss the DWRG method and 
also our strategy. In Section [II, we explicitly perform 
the transformation of the 3D Hamiltonian of Eq. (|l|) 
from the phase to the Coulomb gas representation. Our 
numerical method is explained in Section Finally we 
discuss our numerical results in Section ^ and give a brief 
discussion of some of the effects of weak disorder, a < 1, 
and of finite screening of vortex - vortex interactions. 



II. STRATEGY 

The general idea behind a DWRG is to compute, an- 
alytically or numerically, the energy AE{L) of a domain 
wall in a system of linear size L and fit this to a finite 
size scaling form 



AE{L) - 



(2) 



where is a stiffness exponent, whose sign is of funda- 
mental importance, li 9 < 0, AE{L) vanishes in the 
thermodynamic limit. The energy of the domain wall or 
defect excitation vanishes which implies that, for T > 0, 
the probability of the defect Pl ~ e'^^^^'^l^^ ^ 1 as 
L ^ CO. This in turn implies that the density of such de- 
fects is finite when T > and there will be no resistance 
to an infinitesimal applied force and the system has no 
order. This is analogous to the vanishing of the shear 
modulus in a liquid, the superfluid density in a super- 
fluid or superconductor and the spin stiffness constant in 
an isotropic magnet when T > T^. On the other hand, 
if > 0, such defects will have zero probability when 



L = oo and the system will have finite stiffness and will 
be ordered at sufficiently small T > 0. 

In a uniform system without disorder, the definition 
of the energy of a domain wall of size L, AE{L), is in- 
tuitively obvious. The first step is to find the ground 
state (GS) energy of a system of size L, which requires 
applying boundary conditions (BC) which are compati- 
ble with the GS configuration. For a ferromagnet, this is 
straightforward to implement as the GS configuration is 
known to be one with all spins parallel and periodic BC 
are compatible with this. To impose a spin domain wall 
perpendicular to the x direction, one simply changes the 
BC to antiperiodic along x and periodic in the other d— 1 
directions. Then it immediately follows that 



AE{L) = Eap{L) - Ep{L) ^ L'' 



(3) 



where n = 1 for an Ising model and n — 2 for a system 
with a continuous symmetry such as XY and Heisenberg 
models. 

One would like to use the same strategy for random 
systems described bypEq. (^, as suggested by Anderson 
for Ising spin glasseS. However, it is not so clear how 
to proceed because, for a particular sample (realization 
of disorder), neither the GS configuration nor compatible 
BC is known so computing the defect energy AE{L) is 
problematical. Assuming AE{L) can be calculated, the 
stiffness exponent is defined by the scaling ansatz 



{AE{L))^L' 



(4) 



where (• • •) denotes an average over realizations of disor- 
der. To our knowledge, it is not known how to calculate 
analytically either the GS energy Eq{L) or the energy 
E]j (L) of the system containing a defect relative to this 
GS which means that one must proceed numerically. A 
number of conceptual and technical difficulties are ap- 
parent. The first, and most important, is the technical 
problem of computing the energy difference AE{L) be- 
tween the energies of the system subject to two different 
BC. We ultimately want the disorder averaged defect en- 
ergy {AE{L)) = {Ea(L) - Eb{L)) where Ea{L) is the 
lowest energy of a particular sample subject to BC de- 
noted by a and Et{L) with BC denoted by b. We need 
the individual energies Ea{L) and Eii{L) essentially ex- 
actly because the uncertainty in {AE{L)) must be kept 
as small as possible. Also, to our knowledge, there is 
no proof that the scaling ansatz of Eq. (^ is a correct 
description and, even if it is, the only thing we can be 
sure of is < (d — 2)/2. All results are based on fitting 
data to the scaling form of Eq. (^ so one is attempt- 
ing both to verify the scaling ansatz and to estimate a 
numerical value of 9. For any conclusion to be believ- 
able, the data must have both very small errors and fit 
Eq. (|^) extremely well. The first requirement of highly 
accurate data is the most important as the estimate of 9 
depends on this. Assuming that Ea{L) and Et,{L) can be 
determined exactly for each sample, then AE{L) is also 
known exactly for each sample and the errors in {AE{L)) 
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are 0{N~^^^L'^~^) where A'' is the number of samples of 
size L in d dimensions. If the energy minima Ea^b{L) 
are not found exactly, a crude estimate of the errors in 
/\E{L) is 0{N-^/'^L'^) but this is certainly too low as 
failure of the algorithm to find the true minima because 
of being trapped in a metastable state of energy E > Eq 
will cause systematic errors of unknown magnitude. Em- 
pirically, we find that this can readily cause errors larger 
than {AE{L)) which makes the data point useless. This 
is most likely to happen for large L because the CPU 
time required grows uncontrollably, as do the errors, so 
the large L data becomes unreliable. This technical dif- 
ficulty limits the accessible sizes L to small values as one 
must keep errors in individual data points small. 

We are forced to conclude that the accessible sizes L 
are limited by the necessity of finding essentially exact 
global energy minima of each of a number N of sam- 
ples subject to certain, yet to be defined, BC. To our 
knowledge, there is no algorithm applicable to the sys- 
tems of interest which will find exact minima in polite 
nomial time, such as the branch and cut algorithmo 
for the 2D Ising spin glass or nvmrierically exact combi- 
natorial optimization algorithmsE3 for gauge and vortex 
glass models in the infinite screening limit, so we have 
to live with the fact that our problem is NP complete 
and the required CPU tiHie-|explodes as L increases. We 
use simulated annealingEiH to estimate the lowest ener- 
gies, which seems considerably more efficient than sim- 
ple quenching to T = 0, but we are unable to go beyond 
L = 7 in 3D and L = 10 in 2D. We wish to extract the 
stiffness exponent from the scaling ansatz of Eq. (^) 
with a single power law and this makes sense only if the 
errors on individual data points are very small and the 
fit to the assumed scaling form is extremely good. In our 
opinion, the only sensible strategy is to obtain very accu- 
rate estimates of {AE{L)) for the limited sizes L which 
are feasible for the computer power available. 

In the phase representation of Eq. (]l|), the configura- 
tion space to be searched for the global energy minima 
Ea.b{L) is rather large as the phases 9i G (0, 27r] are con- 
tinuous variables. Searching this space in finite time is 
not feasible as most of the allowed configurations of the 
6i are not even local energy minima. It is well known that 
random XY models with a Hamiltonian of Eq. ([^) can 
be written in a Coulomb gas (CG)nOj| vprtex representa- 
tion via a duality transformationHEZrEa which leaves the 
partition function invariant. This expresses the Hamil- 
tonian in terms of charge or voijtcjtjconfigurations which 
are already local energy minimacjcij. Thus, a reformula- 
tion of the Hamiltonian of Eq. (|l|) as a CG performs a 
partial minimization. A further minimization of the CG 
Hamiltonian corresponds to searching the much smaller 
space of local minima. Reformulating the problem of Eq. 
(|l|) including the BC in CG language is clearly a worth- 
while exercise as it dramatically reduces the number of 
configurations over which we have to minimize, despite 
introducing long-ranged Coulomb interactions between 
vortices. The transformation is carried out in Section III 



for the model of Eq. (|) in 3D. 

The final problem is to define what is meant by a do- 
main wall and the BC needed to induce a wall in a finite 
system of size L vn d dimensions. We imagine the sys- 
tem of Eq. on a torus in 2D or a hypertorus in 3Z?, 
which corresponds to imposing periodic BC in the phases 
^i+Le,, Oi where is a unit vector in the direction 
li^x,y,---,d and i = (i^, • ■ • , i^) with = (1, 2, • ■ • , L). 
The phases at corresponding sites («, j) on opposite faces 
are coupled by some interaction V{9i, 9j, Aij) which may 
be regarded as defining the BC. In principle, the GS is 
obtained by minimizing the energy with respect to the 
L'^ bulk variables 9i and all forms of V . This program is 
beyond our ability and we restrict ourselves to those V 
which induce a spin or chiral defect, which are related to 
the continuous and reflection symmetry respectively. To 
impose a spin defect, we choose V = V{9i — 9j — Aij). 
The plaquettes between the opposite faces are indistin- 
guishable from the others and play no special role. We 
therefore keep fixed the frustrations /r — X]p(r) -^ijl'^'^ 
where the sum is over the bonds in a clockwise direction 
of the elementary plaquette centered at r. We still have 
the freedom to add to every bond in the direc- 
tion between opposite faces which imposes a global twist 
in the phase round a loop circling the hypertorus 
in the direction jj,. This is equivalent to a gauge trans- 
formation Aij — > Aij + L on every bond ij in the 
direction /i. The lowest energy, Eq[\^)^ is 27r periodic 
in A^ with a minimum at some A^ which depends on 
the sample. To induce a spin domain wall normal to x, 
one changes the twists from their best twist (BT) values 
Aj^ — > A^ -f- ttS^^x- The minimum energy subject to this 
constraint gives Esd{L), the energy of the system of size 
L containing an extra spin defect. Note that Ego > Eq 
for every sample but Eq is not necessarily the absolute 
minimum as some other functional form of V may give 
a lower energy. However, even if Eq is not the true GS 
energy but is the energy of a state with some excitation 
from the GS, this method of inducing a spin defect en- 
sures that any excitation in the BT configuration will also 
be present in the state with an extra spin domain wall so 
that AEf'^{L) EE EsoiL) - Eq{L) > is not affected by 
these. It is convenient, but not necessary, to define the 
spin defect energy by a twist of tt from the BT value A|J . 
This choice yields the maximum defect energy AE{L). 
Any other choice < e < tt yields the same spin stiffness 
exponent 9f^ defined by 



{AE^'{L,e))^A{e)L'^ 



(5) 



The size e of the twist from the BT value A|J affects only 
the amplitude A{e) which is a maximum at e = tt. 

A chipJ domain wall is induced by imposing reflec- 
tive BCEj which means that corresponding sites (ij) 
on opposite faces are connected by interactions V ~ 
V{9i + 9j — Aij) which is equivalent to a reflection of 
the spins about some axis. In principle, one follows the 
procedure for a spin domain wall to obtain the chiral 
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defect energy AEc = Ecd — Eq where Ecd is the min- 
imum energy of the system with these modified inter- 
actions V connecting opposite faces. However, there is 
no reason to expect Ecd > Eq as the BC defining Eo 
may trap a chiral defect in some samples and, in such 
cases, the modified interactions V wiU cancel the chiral 
defect to give EcD < Eq. This phenomenon has been 
observed psewiously in numerical simulations of the XY 
spin glassOO. We therefore define the chiral defect en- 
ergy as A£;f'^(i) = \Ecd{L) - Eo{L)\ and the chiral 
stiffness exponent Of^ by a finite size scaling ansatz anal- 
ogous to Eq. (^) 



{AEr{L))-L'' 



(6) 



Note that the Hamiltonian of Eq. (|^) is truly invariant 
under reflection 9i —Oi in the XY spin glass case when 
Aij = 0, TT as Aij = ±7r are equivalent. The Hamiltonian 
with uniform distribution of Aij such as a gauge glass 
is not truly invariant under reflection which would also 
require Aij — » —Aij for the Hamiltonian to be invariant 
as it lacks the reflection symmetry. 

However, in an XY spin glass there are two possible 
types of order, spin glass order and a chiral glass order 
each with their own stiffness exponent of Eq. (|5|) and Eq. 
(^. Recently, an important prediction was made that 
9s = Oc < Q for an XY spin glass in dimension d < di 
where di is the lower critical dimensionO. Although not 
rigorous, the arguments are very plausible and supported 
by analytic calculations on sirnule .one dimensional sys- 
tems in which Og — Oc exactlyoO. There is a notable 
lack of analytic results in this field with which to test nu- 
merical simulations and to our knowledge this is the only 
one existing at present. We have checked our numerical 
method in d — 2 < di and get agreement with the ana- 
lytic prediction that d^-r= 9c — —0.37 ± 0.015 to within 
numerical uncertaintycll Assuming the conjecturaSJ is 
correct, this agreement gives some confidence in our def- 
inition of domain wall energies as discussed above and 
in our numerical method in d = 3 using the CG rep- 
resentation. There is no analogous equality of the stiff- 
ness exponents in a XY spin glass ior d — i > di so 
we do not attempt to estimate Oc in 31? but concentrate 
on the spin stiffness exponent 9s- Also, at present, we 
are unable to derive an expression for a CG Hamiltonian 
with refiective boundary condi|liQi|Lin d = 3. All previous 
workj-ap-t^ XY gauge glassEltlLrlj and on the XY spin 
glassE§E3EZl using the T = DWRG method have used 
different definitions for domain wall energies. Minimiza- 
tion with respect to the global twists is omitted, the 
lowest energy with = is called Ep and the lowest 
energy with A^ = tt is called Eap- Neither of these BC 
is compatible with the GS configuration as both must 
induce some excitation from Eq. Nevertheless, the spin 
defect energy is defined by A£;f^ = \Eap- Ep\ and the 
spin stiffness exponent Of'^ by 



We call this a random twist (RT) measurement as, for a 
particular sample, the twists A^ = 0, tt are two arbitrary 
random choices relative to the best twist A?,, which is 
the twist which yields the lowest energy. In a uniform 
ferromagnet, A° — which is realized by periodic BC 
and A^ + tt by antiperiodic BC. For a particular realiza- 
tion of randomness, A° is the analogue of periodic BC 
in a uniform ferromagnet. 



III. TRANSFORMATION TO COULOMB GAS 
REPRESENTATION 

In this section, we discuss the CG representation of the 
Hamiltonian of Eq. (|l|) including all finite size contribu- 
tions. This representation parameterizes the energy in 
terms of the topological excitations on a torus in 2D and 
a hypertorus in 31? and includes global excitations which 
wind round the whole hypertorus. These latter excita- 
tions are very important for a finite system and are vital 
for finite size scaling considerations when one is limited 
to small system sizes L. Also, every allowed configuratiaa. 
of topological excitations is a local energy minimun£3'c2l 
as spin wave excitations decouple from the vortex excita- 
tions, which allows us to obtain a more accurate estimate 
of energy minima than using the phase representation 
of Eq. (|l|) with the limited CPU time available. The 
transformation of the two dimensional XY model to the 
CG representation including boundary c ontrihu tions has 
been discussed in detail in earlier worksE3E3E3. In this, 
section, we use the method of Ney-Nifle and HilhorstEHl 
to transform to the CG representation in 3-D. 

We first replace the potential V{(j)) in Eq. (|l|) by a 
piecewise parabolic potential which is equivalent to a 
VillainE3 potential at T = 0. The partition function for 
a. L X L X L gauge glass model in 3D is 



Z= r Ylde^Y^ exp -f3J ^ - A 

t {riij} I <ij> 



(8) 



where Oij = 9i — 9j — 2'Knij and where riij — — ".ji 
integers on the bond < ij >. By choosing one phase, 9o, 
as a reference, the partition function can be written as 



~Tiii are 



Z = 



d9Q 



Y[ d9ij cxp 



-/3J(% -Ay)' 



(7) 



^ n ^ ( X! "1°^ 2Tr\ si )9^j mod 2tt j 

^ n ^ ^x '^'^ ^ ^x '^'^ 

Here, r is the coordinate of the center of an elementary 
cube of the original lattice which corresponds to the co- 
ordinate of a dual lattice site, r^y is the coordinate of the 



FIG. 1: The relation between the coordinate of the center of 
cube r (dual lattice site) and the coordinates of the plaquettes, 
Txy, Vyz, and Vzx, associated with the cube at r. i and are 
sites of the original lattice. The random bond variable Aij^ 
is relabelled by and similarly for the others. We assign 
three independent random bond variables, A'^ with fi — x, y, 
z, to each cube at r as shown. 



center of the elementary plaquette in the xy plane and 
similarly for Vyz and Vzx- Note that for a 2D system in 
the xy plane, r^y are the dual lattice sites r. Since each 
cube has six faces (plaquettes), each of which is shared 
by two adjacent cubes, to each dual lattice site r we as- 
sign three independent plaquettes with centers at r^y, 
Tyz, and as shown in Figs. (|l|) and (||). I]p(r,„)% is 
the circulation of 9ij round the plaquette in the xy plane 
of the cube at r, 



jyi 



(10) 



and 9ij is the circulation of 9ij along an arbitrary 
loop in the x-direction round the hypertorus, 



(11) 



where with /i 



lattice site and iy and i 



X, y, z is the coordinate of the original 
are fixed. Other summations 
are defined similarly. Note that one needs to consider 
only one global loop on the hypertorus in each direction. 
Circulations round other global loops can be expressed in 
terms of circulations round any three chosen global loops 
and round elementary plaquettes. It is clear from the 
definition of 9ij and periodic boundary condition imposed 
on the 9i, these circulations are integer multiples of 27r. 



FIG. 2: Graphical explanation of our symbols, r denotes the 
center of the cube and light solid lines join original lattice 
sites. 



Since the delta functions can be rewritten as follows, 



S j 9ij mod 27r 



1 

2^ 



oo 

— — oo 

oo 



P{^xy, 



S I 9ij mod 2tt \ — — exp 
the partition function now becomes 

-/?J(% 



/oo 
n exp 



X exp 



X exp 



E(< E 

Piry.: 



P{r.y) 



9,j + iuy 9ij +inzJ2^ 



(12) 



where n!^ and with /j, = x,y,z are integers, and X)n = 
and X;{„,} = E„!^ The sum Y.r is over the 
dual lattice sites r at the centers of the elementary cubes 
and is an unimportant constant. We use the notation 
Sp(r ) ^o denote a sum over the bonds, in a clockwise 
direction, of the plaquette in the xy plane centered at 
Vxy and to denote a sum over bonds on a global loop 
round the hypertorus in the x direction. To perform the 
integrations over {%}, we choose the three global loops 
around the hypertorus as shown in Fig. (^. To deal 
with these global loops mathematically, we introduce the 
following quantities, 



= 



1 if r = (a;, 1, 1) with x 
otherwise 



1,2,- 



.L 



(13) 
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similarly for and 6^. For example, the cubes at r = 
(x, 1, 1) have a part of the global loop in the x direction. 

The definition of Af^ associated with the cube at r is 
also shown in Fig. (|l|). As the plaquettes, one can assign 
three independent to each cube. After relabeling Aij 
by A^, performing the integrations in Eq. ( |l2|) over 9ij, 
the partition function becomes 



Z = Zq cxp 



4/3J 



X exp 



(14) 



We use the following notation for the discrete derivative, 

(V X nr). EE [nl ~ nl^^) - (n^ - nl_^) (15) 

and similarly for other components of V x rir. Note that, 
when r — (1, y, z), then r — x = [L, y, z) due to periodic 
boundary conditions, and similarly for r — y and r — z. 
Applying the Poisson summation formula to 

E E / due'^^'Viu). (16) 

n— — oo q— — OQ 

the partition function of Eq. ([T^ ) becomes 

/oc 
JJ^ dur exp 2TTi qr • Ur 
II .j^q^^ oo J. J. 



X exp 



X exp 



^E E {(VxUr)„+<S>J^ 



4/3J 



r a.—x,y.z 



-*E E {(VxUr)„+<5?"a}^? 



r a—x,y,z 



(17) 



It is convenient to change integration variables to to 
make Eq. ( p7| ) more symmetric 



(18) 



where r = {x,y,z). The partition hmction of Eq. (|l 
now becomes 







J/ 






= — 




X 




= Ur - 


X 

— n,, - 
2i ^ 


2/ 



2L 



/OO 
n^VrCXp 

^^qr-(rxn} 



X exp 



X exp 



X exp 



4/3J 



r a— 



^E E {(VxVr)„ + !^}A^ 



r a— a;,y,2: 



(19) 



The terms linear in in the exponent of Eq. 
( |l9| ) vanish due to the periodic boundary condition 

J2r (^r ~ "r-y) ~ 0. Eq. (|l9| ) cau bc simplified by 
introducing frustration variables at site r 



(20) 



from which and are obtained by cyclic permutation 
of xyz. After some algebra, the partition function of Eq. 
( p^ ) becomes a form suitable for integration over 



Z 



27ri 

2r 



|qr ■ (r X n) I 



X exp 



X exp 



4/3J 



r a.—x,y.z 



27ri^(qr - fr) • Vr 



(21) 



To evaluate the integrals over u{f in Eq. (pi| ) it is con- 
venient to take the Fourier transform 



Vr = L-i^e^''"'v(k) (22) 



where k = {kx,ky,kz) and ki — •^Trii with rrii — 
0, 1,---,L — 1. We also decompose v(k) into longitu- 
dinal, Vi(k), and transverse, vx(k), components where 



(1 



)(1 



^2(k) = E 



■v^ik) 



Aw = 



6 — 2 cos k^ — 2 cos ky ~ 2 cos 



(23) 



define the longitudinal and transverse projection opera- 
tors Lapik) and TQ^(k). In Fourier space, the integration 
over Vr in Eq. ( pT] ) is straightforward 



J|£ivi(k)dvT(k) exp 



X exp 



27ri ^ { PT (k) • vj, (k) + Pi (k) • Vi (k) } 



(24) 



where p(k) = q(k) — f(k). Integration over VL(k) gives 
PL(k) — 0. In real space, this is the condition that the 
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discrete divergence of the vorticity (charge) qr at any 
dual lattice site r obeys V • Qr = since V • fr = 0. 
The integration over the transverse components VT(k) 
are simple Gaussian integrals and are easily performed. 
Integration over Vi(0) and vt(0) yield the neutrality 
condition 



p(0)=^Pr 







(25) 



The final step in this rather technical derivation of the 
Hamiltonian in the CG representation is to apply the 
Poisson summation formula of Eq. (n&j to eliminate the 



in favor of global vortices or charges q^i . The Gaus- 



sian integrations yield the partition function 



Z 



f3H{qr,^v,q^aJ^.l)] 



(26) 



The Hamiltonian H is identified as 



H = (27r)V^(qr-fr)-(qr' -fr')G(r-r') 



J r TT 
2ZIZ 



J r TT 
2LII 



yPr 



XPr-ZPr) +Qv] 
r 



(27) 



where G{r) ~ L^^ Sk5^o(^^P(*^ ■ r) — 1)/Ak is the lattice 
Green's function on a simple cubic lattice in 31? with 
periodic BC and the quantities are 

= 7r^(zp^5y4 - ypl5,,i) + 2^i(g,i - U) (28) 



with Qy and Qz obtained from Qx by cyclic permutations 
of xyz. In Eq. (p8|), fxi is the circulation of Aij along 
the chosen global loop round the hypertorus in the x 
direction 



27r/:rl = X!^r=(x,l,l) 



(29) 



and similarly for fyi and fzi- The integers g^i are in- 
terpreted as circulations of the phase round the three 
independent global loops encircling the hypertorus. Eqs. 
(|^) and (11) are the main results of this section. These 
expressions give the energy of a system of size L on a 
hypertorus. To find the minimum energy Eq the Hamil- 
tonian H = H{qy., fr, (7^1, /pi) is minimized with respect 
to the bulk integer valued vector charges qr, the inte- 
ger valued global winding numbers g^i and the global 
frustrations /^i. In the case of the XY spin glass the 
9^1 ~ /^i of Elq. ( p8| ) are restricted to be integer or half 
integer while for the gauge glass this can have any real 
value. Note that minimizing with respect to q^i and fxi 
is exactly minimizing with respect to the twist A^; and 




FIG. 3: The L x L x L system is represented by a cube. The 
thick lines are our choices of the three global loops around 
the whole system. 



the best twist corresponds to /°i, the value of /^i 
at the energy minimum. Adding a global twist A^ to 
the phases is exactly equivalent to changing the global 
frustration from its original value /^i —>■ ff^i -f Ap/27r. 
As discussed in Section ||, a spin domain wall is induced 
by A° — + A|^ -f TT which, in the CG representation, is 
fxi fxi + 1/2- Note that the frustrations fr are given 
in terms of the random Aij and are kept fixed during the 
minimization. 

To confront our nupierical results with the only existing 
analytic predictiortJ we need expressions for both spin 
and chiral domain wall energies in dimension d = 2 < di. 
The CG representation of the XY spin and gauge glasses 
including all finita,size corrections have been known for 
some time in 2Z3E3Ea and, for the sake of completeness, 
we quote the necessary results below. To study the spin 
stiffness, we join corresponding sites on opposite faces by 
V — V{di — 9j — Aij) and the CG Hamiltonian is 

H = (27r)2j^(qr-/r)G(r-r')(gr' -/rO 

rr' 

J 



(30) 



where 



= -277 



= -27T 



L{qxi 
L{qyi 



/xl) +^(gr - fr)y 
r 



G(r) 



-y 



1 



4 — 2 cos kx — 2 cos k,. 



(31) 



Here, as in 31?, r — {x, y) represents the coordinates of 
the dual lattice sites and G(r) is the lattice Green's func- 
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tion with periodic BC on the 9i. From Eq. (|8|), we see 
that the difference of a factor 2 in the prefactor of Eq. 
( ^ ) from other works is in the couphng constant J. Also 
note that the Hamiltonian of Eq. ( ^o|) describes the XY 
spin glass when the frustrations /r and the global frus- 
trations f^i are restricted to (0, 1/2). In the gauge glass, 
they can have any real value. As in 3-D, the topological 
charges are integers as are q^i because of the periodic 
BC in the 9^. 

The last piece of information we need is the CG repre- 
sentation of the Hamiltonian of the XY spin glass in 2D 
with a chiral domain wall imposed. Tiiis is to be found 
in the paper of Ney-Nifle and HilhorstL3. A single chiral 
domain wall is induced by joining opposite faces by in- 
teractions V = V{9i -1r-ftn— ^ij) which is equivalent to 
imposing reflective BCx3Ej. In turn, this is equivalent to 
doubling the size in (say) the x direction to a 2L x i lat- 
tice in which the extra half is a charge conjugated image 
of the other. This-|System has two chiral domain walls 
with HamiltonianEJ 

Hr = 27rV^(<Zr - /r)G(r - r')(gr' - A') (32) 

r,r' 

where G(r) is the Green's function for a 2L x L square 
lattice with periodic BC and also with qr+Lx = ^qr and 
fr+Lx = —fr- Note that the sign reversal of the frus- 
trations /r is not necessary for the spin glass because 
fr = ±1/2 are equivalent. 

The form of the energy of Eq. (|3|) is used in the sim- 
ulations to estimate 9c the chiral stiffness exponent as it 
is intuitively more transparent and more convenient than 
the corresponding expression with a single chiral domain 
wallEj. Unfortunately, we have been unable to derive the 
analogous expression in 311. to Eq. (|3^) so we have no 
independent estimate of 9jE3 for the 3D XY spin glass. 



IV. NUMERICAL METHOD 
A. Minimization Algorithm 

In Section || we argued that it is necessary to find the 
energies Eo{L) and Ed{L) essentially exactly for every 
sample to control the errors in the small domain wall 
energy {AE{L)). To our knowledge, for the systems of 
interest no algorithm exists which can locate the global 
energy minima in polynomial time. We are left with two 
methods: (i) repeated simple quenches from an arbitrary 
initial configuration to T = followed by a downward 
slide to the nearest local minimum and (ii) simulated 
annealingo which is considerably more efficientt3. By 
this we mean that, for the same CPU time, simulated 
annealing finds a lower energy than simple quenching. 

We start with the system in some randomly chosen 
configuration and quench to some Tq, determined by trial 
and error, for each size L. Then we do a Monte Carlo 
(MC) sweep through the system by inserting charges 



q = 1 and q' = —1 on an arbitrary pair of nearest 
neighbor sites of the dual lattice in 2D and accepting or 
rejecting the move according to conventional MC rules. 
This has the effect of inserting new charges, annihilat- 
ing charges or moving a charge by one lattice spacing 
while maintaining charge neutrality 9r = 0. In 31?, 
the elementary excitation is a loop of charge round an el- 
ementary square with vertices at dual lattice sites. This 
maintains V • qr = 0. The closed loop of charge can lie in 
any of the three orthogonal planes of the cubic lattice and 
the charge can circulate round the loop in either direc- 
tion, making 6 possibilities with the center of the loop at 
a fixed but randomly selected position. The temperature 
T is then reduced to Ti = oTq with a < 1 whose value 
is again determined by trial and error and the procedure 
iterated a large number, N, of times to reach a lowest 
temperature Tn = a^Tg « 0. Of course, the system 
may be trapped in a deep metastable well with barriers 
too high for the MC passes to overcome so the whole 
annealing sequence is repeated M times from different 
random initial configurations and the lowest energy out 
of all the NM trials recorded. Again, this does not guar- 
antee that the global minimum energy is found but this 
method does have a few checks built in. At the crudest 
level, the best twist condition Esd{L) > Eo{L) must be 
obeyed for each sample since Eo{L) is, by construction, 
the lowest energy of the system subject to the BC given 
by the interactions V — V{9i — 9j — Aij) across oppo- 
site faces. Esd{L) is obtained by /"^ f^^ + 1/2 where 
f^i is the value of /ui which makes the boundary terms 
of Eqs. ( p7| ) and ( |30| ) vanish, corresponding to the best 
twist A". It is clear that Esd{L) is the energy of the 
system containing an extra spin domain wall compared 
to the system with energy Eo{L). As discussed earlier, 
Eo{L) is not necessarily the true GS energy as the sys- 
tem may contain some chiral domain walls. However, by 
construction, Esd{L) is the energy of the system with 
the same chiral defects and an extra spin defect. 

If any sample violates the BT condition, clearly the an- 
nealing is not sufficient and one can either increase the 
number MN of annealing attempts or just discard that 
sample. Increasing MN involves a significant increase in 
CPU time particularly for the larger sizes L so the choice 
depends on the time available. However, even if the BT 
condition is satisfied for every sample, there is no guaran- 
tee that the lowest energies found are true global minima. 
To improve the chances that true minima are achieved, 
where possible we did two independent simulations on 
identical samples with different pseudo random number 
sequences and, if the same energy minima are found in 
both simulations, this is defined to be the true minimum. 
This procedure is very expensive in CPU time and our re- 
sources did not permit this last check to be performed for 
every sample, particularly for the largest systems. This 
check was done for at least a few randomly selected sam- 
ples, except for i = 7 in the 3D gauge glass. Averaging 
over disorder was performed over as many samples as 
possible with the aim of making the uncertainty in the 
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mean domain wall energy {AE{L)) less than 3% which 
requires averaging over at least 10'^ samples. The first 
set of checks are done on every sample which makes it 
fairly probable that exact minima Eq{L) and Eo{L) are 
obtained. We can assume with some confidence that the 
uncertainty in {AE{L)) is purely statistical and 0{N^^'^). 
Averaging over about 10"^ samples leads to an acceptable 
uncertainty of about 3% in {AE{L)). 



B. XY Spin Glass in 2D 

This particular system with a Hamiltonian of Eq. (|^) 
in 2D is not particularly interesting in the sense that 
there exists no finite temperature transition. It has been 
known for a long time that di > 2 and both spin and 
chiral stiffness exponents are negative. However, the 
situation has been controversial due to the possible ex- 
istence of different stiffness or correlation length expo- 
nents for spin and chiral glass order in 2D. Previous 
estimates of the spin and chiral stiffness exponents are 
summarized by 9s ~ 29c ~ —0.78 based on extensive 
numerical simulations piehpas_D|WRG and finite temper- 
ature MC simulationsEaOOBEl. However, Ney-Nifle 
and HilhorstEy made a non-rigorous but very plausible 
conjecture based on analytic considerations that for di- 
mension d < di, 9s = 9c < 0. This is supported by exact 
analytic resultS|JMi simple models (i) a XY spin glass on 
a ladder latticeEa where the common correlation length 
exponent — \9\^ — 0.5263... and (ii) the XY spin glass 
on a tube lattices with v = 0.5564.... The key disagree- 
ment between numerical estimates and analytic theory is 
in the conjectured equality of 9s and 9c- If one makes 
the hypothesis that earlier estimates are in error, then 
the most likely reason is that the simulations are com- 
puting the wrong quantity when estimating one or both 
of the stiffness exponents. It is extremely unlikely that 
all the simulations are in error for simple technical rea- 
sons as all find the same values of 9s and 9c, but 9s 9c- 

To test this hypothesis, we have performed simulations 
of the 2DXY spin glass in the CG representation us- 
ing Eq. ( pO[ ) to estimate the spin stiffness exponent 9s 
using both BT and RT measurements with the results 
9f^ = -0.37 ± 0.015 and 9^^ = -0.76 ± 0.015EI. These 
numbers were obtained from sizes L — 4,5, 6, 7, 8, 10 av- 
eraging over 2560 samples for i < 8 and 1152 samples 
for L = 10. As expected, the value of 9f^^ agrees with 
all previous estimatesc^O, all of which use the RT mea- 
surement in some form. Both the BT and the RT data 
fit the scaling ansatz of Eq. (jj) equally well and some 
other information is needed to decide which value of 9s , 
if cither, is correct. Both cannot be correct as both are 
supposed to measure the same quantity. The necessary 
information is in the chiral stiffness exponent 9c which 
we measure by simulating Eq. (^) on a 2i x i sys- 
tem. Again, both BT and RT measurements were made 
with the same range of L and the same number of sam- 
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FIG. 4: Top to bottom; L dependence of AE^"^ , AEf^ , 
A_Bf ^, ASf ^ for the 2DXY spin glass. 



pies as for 9s with the-J-esults 6lf ^ = -0.37 ± 0.010 and 
6»f^ = -0.37±0.015.EI At first sight it is surprising that 
both measurements give the same value for 9c to within 
numerical uncertainty while the values of 9f^ and 9^^ 
differ by a factor of 2. Note that the boundary terms a a 
of Eq. ( pO| ) which contain the twist parameter /^i van- 
ish in the BT condition and such boundary terms are 
absent from Eq. (32). Thus, any measurement with 
a CG representation in which boundary contributions 
are absent is automatically a BT measurement. Since 
9c = 9f^ « -0.37 and 9c ^ 9^^ « -0.76, assuming that 
the conjecture is correct, we conclude that the BT mea- 
surement yields a reasonable estimate of the true 9s while 
the commonly used RT measurement yielding 9s ^ ^ 9c 
is not an appropriate method for the small values of L 
accessible at present. Our simulation results are shown in 
Fig. (|^). We do not understand what, if anything, 9^^ 
means despite the apparent excellent fit of {AEf''^{L)) 
to the scaling ansatz as this is the energy difference of 
the system subject to two random choices of BC and we 
can see no reason why it should scale as L^' over decades 
of L. 



V. RESULTS 

Our investigation of the 2DXY spin glass establishes 
that the BT method of measuring domain wall ener- 
gies can rearoduce the one and only available analytic 
predictionE3 which is relevant for our purposes. This 
gives some encouragement to venture into areas where 
no such analytic guide exists. In this Section, we report 
some new results on the spin stiffness exponent 9s for the 
XY spin gla ss in iD (V A), and the gauge glass in both 
2D and 3D ( VB ). We also perform simulations on sys- 
tems with varying strengths of disorder ( |V C| ) where we 
study the renormalization group flows for both the cou- 
pling constants J{L) and disorder strength A{L) with 
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increasing length scale L. All our results are numerical 
and the exponent 9s = Of^ is estimated by fitting esti- 
mates of {/S.Ef^{L)) to the scaling ansatz of Eq. (^). 
The system sizes L are very small as they are severely 
constrained by the necessity of controlling the errors in 
both Eo{L) and Ed{L), which are estimated by indepen- 
dent simulations. However, the fit to the scaling ansatz 
is very good despite the very few data points for all cases 
we have studied and is comparable with the same pro- 
cedure carried out on the 2D Ising ferromagnet. In this 
case, simulations for sizes L = 2,3,4,5 are sufficient to 
reproduce the exact = 1 to high accuracy. 



A. XY spin glass in 3D 




This system has been somewhat controversial for some 
years and is not yet settled. It has been believed that 
that di > A for spin glass order, and d; < 3 for chiral 
glass order and that, in 2D and 3D, AomAndL/shitfl vari- 
ables decouple and order separatelyuEaESEZleS'CJ. This 
allows for the widely accepted scenario that in SD spin 
glass order sets in at Tsg = whereas chiral glass or- 
der sets in at Tcg > fhrThis scenario is based on MC 
simulatipiis at finite Tc3c3 and on T = defect energy 
scalingEj'cEl using the RT metiiod. Attempts to show rig- 
orously that /Tfi i g i t in 31X5 fail if reflection symme- 
try is brokenr ttttI at finite T. The first cracks in this 
widely accepted sceBario appeared recently when Mau- 
court and GrempelEZi published the results of a large 
scale defect energy scaling study of the 3D XY spin 
glass model of Eq. ([^). They used the RT method, 
as all previous defect energy scaling studies have done, 
with sizes L < 12 in 2D and L < 8 in 3D. Although 
the fit to the scaling ansatz is not good due to strong 
crossover effects and large uncertainties in the large L 
data which was used to estimate 0^^,it is clear that both 
(A£;f^(i))and {AE^'^{L)) are increasing with L which 
implies that there is both spin and chiral glass order at 
sufficiently small T > 0. However, as argued above a 
valid numerical method must|4j"ield 9^ = 9c in 2D while 
they obtaia_6L « 29c ~ — 0.78E3 in agreement with other 
estimatescJO. Finite T simulations seem to suffer from 
severe equilibration difficultiesEj'LJ which makes any con- 
clusions from them also suspect. 

In view of the lack of reliable results for the 3DXY spin 
glass, we have done some preliminary simulations on very 
small systems with L ~ 2,3,4,5,6 in the CG represen- 
tation with Eq. (H) where /^i = 0, 1/2 in E q. (H). 
Following the method outlined in Section [VB| we esti- 
mate 9f^ = +0.10 ± 0.03 while the data for {AE^^{L)) 
is clearly decreasing with L for L — 2, 3_4j with roughly 
the same slope as found by KawamurECJ for the same 
range of L. Our value 9s ~ 0.10 ± 0.03Eil is to be com- 
pared witlijhe recent estimate 6'f — -1-0.052 ± 0.03 for 
L = 5, 6, TEZI. Although these values are close numerically 
and are equal to within error bars, we are more inclined 
to believe in the former number as we consider 9^^ to 



FIG. 5: L dependence of AE^f^ and A_Ef ^ in 3D. The error 
in the L = 6 point is due to rather few samples. Solid line is 
power law fit and the dotted line is a guide to the eye. 



give the true spin stiffness exponent and we speculate 
that further points will also lie on the scaling ansatz and 
will significantly reduce the 30% uncertainty. The results 
are shown in Fig. (^). Unfortunately, we have no esti- 
mate for the chiral stiffness exponent 9^^ for the 3DXY 
spin glass as we have been unable to derive the 31? ana- 
logue of Eq. dH) to estimate AE^'^{L) and AE^^{L). 
However, because the numerical values of 9f^ and 9^^ 
are the same in 2D and equal to 9s as required, we-^ifa 
reasonably safe in assuming that existing estimateaHj'EII 
of 9c in 3-D are fairly accurate. Even though our spin 
stiffness exponent 9f^ « -1-0.1 has a rather large uncer- 
tainty, it suggests that the lower critical dimension for 
spin glass order is < 3, as for chiraL order. This is to 
be expected from analytic argumentsESI. 



B. Gauge Glass in 2D and 3D 

We have also performed simulations on the gauge glass 
in the CG representation using Eq. (|3^) in 2D and 
Eq. ( p7| ) in 31?. The only differences to the spin 
glass are in the values of f^i and /r which can have 
any value in the interval [—1/2,1/2). The frustrations 
/r — — X)p(r) ^ii/^'"^ correlated random variables 
as the Aij are the independent random variables. Simi- 
larly, in 3D, the frustrations fr of Eq. ( pO| ) are correlated 
random variables with each component in the interval 
[—1/2,1/2). The three global frustrations f^i can take 
any value in the same interval. In both spin and gauge 
glasses, the vorticities qj^ and g^i are integers. We have 
not computed 9c in either 2D or 3D for the gauge glass, 
mainly because we have been unable to obtain the ap- 
propriate expression for the Hamiltonian in the CG rep- 
resentation with reflective BC in 3D and we are unable 
to understand what information such a simulation would 
yield. The major reason for doing such a simulation in 
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FIG. 6: Size L dependence of domain wall energy in 2D. Both 
RT and BT measurements are shown. Solid hnes are power 
law fits. Error bars are not shown if smaller than symbol size. 



the 2D spin glass case was to confirm that our procedure 
is a useful numerical method to estimate the actual spin 
stiffness exponent 6s for both spin and garua&iglasses in 
2D and iD. To compare with earlier workBQBo we have 
also estimated by keeping /^i fixed at some random 
value during the minimization, exacthi, as for the spin 
glass. The results of the simulationscS in 2D for sys- 
tem sizes L — 2,3, 4, 5, 6, 8, 10 are shown in Fig. (||) for 
{AE^'^{L)) and (AE^'^ (L)). Averages were performed 
over about 10^ samples for ea ch siz e L. In this case, all 
the checks discussed in section IV A were performed so we 



assume the errors are purely statistical. Both fit well to 
the scaling ansatz of Eq. with very similar errors with 
the values Of^ = -0.36±0.010 and Of^ = -0.45±0.015. 
ThjELjlatter value is consistent with all earlier estimates of 
^sHfla which is not a surprise as all these were done us- 
ing the RT m ethod in some form. However, as argued in 
section [V B , this is not an accurate estimate of 6^ so this 
number is not the T — spin stiffness exponent. A more 
accurate estimate of this is Of^ — —0.36 ± 0.010 which is 
significantly larger than df^^ , so that the 2D gauge glass 
has a longer correlation length ^(T) ~ T^-i/l^^l than pre- 
viously thought. 

We have also obtained, some estimates of 9s in 3D by 
performing simulationscS of the gauge glass in the CG 
representation using Eq. ( p7| ) with the distribution of 
frustrations fr appropriate for this case as determined by 
taking the Aij uniformly distributed in (— tt, -I-tt]. The 
system sizes are 2 < L < 7 with disorder averaging over 
10^ samples for L < 5, 300 for L — 6 and 60 for the 
largest size L — 7. The uncertainty in {AEf^^L)) for 
L = 7 is very large, but this data point is included to 
check that it is consistent with the behavior deduced from 
the more reliable data of the smaller sizes. The results 



are shown in Fig. (|^) for {AEs{L)) for the unscreened 
gauge glass using both BT and RT measurements and for 
a gauge glass in 3D with screened interactions, A < oo. 
The BT data fit the scaling form of Eq. (^) very well 
for sizes L < 6 with exponent Of^ = +0.31 ± 0.010. If 
the L — 7 point is also included in the fit, we obtain 
0^^ = -1-0.30 ±0.015. These errors in Os come from 
a naive least squares fit of the data to a straight line 
and should not be taken too seriously. The L = 7 data 
is suspect because 4 samples violated the BT condition 
AEf'^{L) > out of a total of only 64, despite running 
highly vectorised code for a few thousand CPU hours on 
a Cray J90. Time did not permit any further checks for 
attaining the global energy minimum for L = 7. We 
cannot be sure that the remaining 60 samples which did 
not violate the BT condition reached their global energy 
minima nor that the energies Esd are determined suffi- 
ciently accurately. The error bar on the L = 7 point in 
Fig. (0) assumes that the uncertainty in {AEf"^{L = 7)) 
is purely statistical and the true uncertainty is proba- 
bly much larger. At least an order of magnitude more 
CPU time is needed for sufficient annealing to reach the 
true minima and to perform the additional simulations 
with different random number sequences to check that 
the minimization algorithm is successful. This is just for 
a single batch of 64 samples and to reduce the uncer- 
tainty to 3%, yet another order of magnitude of CPU 
time would be needed to average over 10"^ samples. This 
is totally out of reach with the computing resources avail- 
able to us. What data we have is entirely consistent with 
the scaling form of Eq. with ds ~ +0.30 with no sign 
of any deviation from this form. 

The behavior of {AEf'^{L)) is also shown in Fig. (|^) 
for sizes L < 6 which is very much like the data obtained 
by earlier simulations. This clearly does not fit the scal- 
ing form of Eq. (^) for these small values of L, but if one 
insists on extracting a value of 6^^ from tJjie data, one 
obtains consistency with previous estimategS for the spin 
stiffness exponent 6*^^ « -1-0.05 ± 0.05. As can be seen 
from Fig. (^), this estimate has no meaning as the data 
clearly does not obey the scalins, ansatz. In fact, as no- 
ticed by Maucourt and GrempeB, (AEf^{L)) seems to 
start increasing with L for L > 5 but, as we argue earlier, 
this may, or may not, eventually scale as {AEf'^{L)) for 
sufficiently large L. Speculation along these lines is fruit- 
less until computers which are many orders of magnitude 
faster become available or until an analytic solution is 
found. We conclude that, for the unscreened gauge glass 
in 3D, the spin stiffness exponent 9s = -t-0.31 ± 0.010, 
which is considerably larger than earlier estimates and 
indicates that the lower critical dimeijaiaia d; < 3. This 
is consistent with finite T MC resultaJ'ElEj for the gauge 
glass in 3D which indicate that Tc = 0{J). This value of 
Tc is very difficult to reconcile withj-jtjhjepSuggestion that 
di ~3 from previous DWRG studieaj'BEifl as this implies 
a very small value of T^j J . 

We have also studied the effects of screened vortex- 
vortex interactions on the spin domain wall energy using 
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FIG. 7: L dependence of domain wall energy in 3-D. Bot- 
tom curve is RT measurement for unscreened interaction. All 
others are BT measurements. Topmost curve is unscreened 
case L — 2 ~ 7. Other curves are screened interactions with 
A decreasing from top to bottom. Solid line is a power law fit 
and dotted lines are guides for the eye. 



the BT measurement in 31?. Screening of the Coulomb 
interaction of vortices is implemented by adding a niass 
term to the denominator of the Green's functioned 



^3 g _ 2cosfca; — 2cosA;„ — 2cosfc2 + A ^ 

_ (33) 
The results are also shown in Fig. (j^). We average over 
10'^ samples for L = 2,3,4 and 250 for L — 5 for several 
values of the screening length A. Screening is clearly a 
relevant perturbation when A is finite and < but 
our small sizes do not permit an estimate of the value of 
9f^- For large screening lengths, {/S.Ef'^{L)) seems to 
scale the same as for the unscreened case but we expect 
that {AEf'^{L)) will decrease as with < at 
length scales which are beyond our computing power for 
any A < oo. These results are consistent with those of 
Bokil and Younglj who studied the question of screening 
using a RT measurement and with Kisker and RiegerEj 
for very strong screening. 



C. Varying disorder strength in 2D and 3D 

We have also performed simulations with various 
strengths of disorder in the CG representation using Eq. 
@ for the 2D case and Eq. (|2^) for 3D where the 
random bond variables Aij axe independently uniformly 
distributed in the range [—air, air) with < a < 1 so 
that (Aij) = and = a7r/2. Physical realiza- 



FIG. 8: KG fiow in 2D. The flows are from right to left for 
a > 0.4 and from left to right for a < 0.35. 



tions of this model are, e.g., an ACY-piagnet with random 
Dzyaloshinski-Moriya interactiong^^Laijd| Josephson junc- 
tion arrays with positional disorderOta where both the 
effective coupling constant J{L) and the effective disor- 
der strengtlL A(-L) at length scale L play a role. Studies 
in 2ii3'EZ[E3 suggest that weak disorder (a > 0) does 
not affect the existence of an ordered phase at interme- 
diate temperature but there is a re-entrant transition to 
a disordpjedfpliase at low temperature. However, recent 
analyticE£(E3tj'E3 and numericalcl studies show there is 
an ordered phase for T < Tc{a) with Tc{a) > for 
< a < Qfc- To study how the disorder strength be- 
haves as the length scale L varies, one must, if possi- 
ble, identify the scaled disorder strength A{L) with some 
measurable quantity. An identifi«Ltion has been pro- 
posed in a recent numerical studycl where the effective 
disorder strength is defined as 2ttA{L) = (|A*'(L)|) with 
A{1) = a/A so that one can follow the flows of both J{L) 
and A{L) with increasing length scale L. With this defi- 
nition, < A{L) < 1/4. A°(L) is the global phase twist 
minimizing the energy of a system of size L for a partic- 
ular realization of disorder. For two phases with energy 
Ei2 = V{6i — 92 — A12), the minimum is at Ox — 62 = ^12 
which is satisfied by applying a "global" phase twist of 
Ai2. Hence, follows the definition of A{L) as a measure 
of disorder at scale L. r-, 

Since the numerical studyO used the Hamiltonian of 
Eq. in the phase representation, we re-investigate this 
model in the CG representation. The simulations were 
performed for L ^ 2,3, 4, 5, 6, 8, 10 in 2D, i = 2, 3, 4, 5 
in 3D and averaged over at least 10'^ samples. The re- 
sults are shown in Fig. (^) in 2D and Fig. (||) in 
3D. We are interested in the stable fixed point values 
at L — > 00 J* and A* as these determine the nature of 
the phases. In 2D, weak disorder {a < ac~ 0.37) seems 
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FIG. 9: RG flow in 3D. The flows are from left to right for 
all a. 



to be marginal and the system seems to iterate toward 
a glass phase with quasi-long range order characterized 
by (J*, A*) where the fixed point values J*, A* are finite 
and depend-Qn .a .. This is consistent with recent ana- 
lytic studicsoElLdcj. On the other hand, systems with 
strong disorder, a > ac^ seem to flow to a disordered 
fixed point (J*, A*) = (0,1/4) which corresponds to a 
non-superconducting glass. 

In 3-D, for a < ac~ 0.57, the system flows to a strong 
coupling limit J* — oo. The disorder strength A{L) ap- 
pears to flow to a finite fixed point value A* which de- 
pends on a. However this is not conclusive from our 
simulations as only very small values L = 2, 3,4, 5, 6 are 
used. This can be interpreted as the zero field version 
of a Bragg glass phase. For large disorder, a > ac^ 
{J{L), A{L)) seem to iterate to their maximum values 
of {oo, 1/4) corresponding to the gauge glass fixed point. 
It seems that, in the absence of screening, A — > cxd, there 
are two glassy superconducting phases at T = which, m 
an applied magnetic field, correspond to a Bragg glassEll 
for a < ac and to a vortex glassEj for a > Uc- 



VI. DISCUSSION AND SUMMARY 

In this paper, we re-investigate the possibility of an 
ordered phase at small but finite temperature T by a nu- 
merical domain wall renormalization group method in a 
disordered XY model in 2D and 3D described by the 
Hamiltonian of Eq. ([|) in the Coulomb gas represen- 
tation. For the ±J XY spin glass in 3D, our simula- 
tions yield the spin glass stiffness exponent6'f'^ w -1-0.10 
which suggests its lower critical dimension is d/ < 3. This 
value of is very different from existing estimates pi 
the chiral glass stiffness exponent in 3D 6c ~ -I-0.47E3 



and Be - +0.56 ± 0.180. The difference between 9^ 
and 9c seems to support the decoupling of two degrees 
of freedom in 31?. For the gauge glass, we estimate 
the stiffness exponent 9^^ = —0.36 ± 0.01 in 2D and 
9f^ = -1-0.31 ± 0.01 in 3D, which are considerably larger 
than all earlier estimates. The latter value is consis- 
tent with^c/ J ~ C(l) from finite temperature MC 
studiesoO^Ej and also strongly suggests di < 3. The re- 
sults for the XY spin glass in 3D are consistent with spin 
glass order at-|21,|>, which iSrin contradiction with all 
other studicscJcJcj except oncE]. 

We also studied the effects of varying the disorder 
strength. In 2Dp-<wij--aHHulations imply that weak dis- 
order is margina]ll3E3E2l'EHl and a system with strong dis- 
order flows to a disordered fixed point. There is no sign 
of a re-entrant transition in our simulations. In 3D, weak 
disorder has little effect and the system flows to an or- 
deredjjhase which is the zero field analogue of a Bragg 
glassEiJ. For strong disorder, the system seems to flow 
to a gauge glass fixed point. The disagreement between 
the stiffness exponent 9f^ and previous estimates is be- 
cause these measure 9f-^ whose meaning is less clear. The 
quantity AD^'^(L) seems more likely to suffer from large 
corrections to scaling as seen in Fig. (^, especially for 
the small system sizes L which are possible to simulate 
at present. However, we conjecture that both measure- 
ments would coincide if much larger values of L could be 
reached. Since our simulations are also limited to very 
small sizes L, it is not possible to draw any definite con- 
clusions from them and more studies are needed to settle 
these problems in random systems more satisfactorily. 

One interesting conclusion we can reach concerns the 
Bragg and vortex glass states in disordered supercon- 
ductors in an applied magnetic field. Recently, one of 
usEj studied the model of Eq. in the strong screen- 
ing limit A ^ in 3D and found that two phases exist 
at T = in the presence of an applied external field. 
The low field, small disorder phase has a well ordered 
vortex line lattice as a ground state with stiffness expo- 
nent 9s = -1-1.0, whereas the high field, large disorder 
ground state is a disordered entangled vortex configura- 
tion with 9s ~ — 1.0l3Ej. We identify the low field state 
as a Bragg glass and the high field state as a disordered 
entangled vortex liquid. In this limit, the evidence is 
strongly in favor of a direct, disorder or field driven tran- 
sition from a superconducting Bragg glass to a normal 
non-superconducting phase. This scenario seems to be 
favored by recent experiments. 

In the absence of screening of the vortex-vortex inter- 
actions, the picture which results from this work is some- 
what different, although the studies here are all done in 
zero applied field. One may argue that increasing the 
disorder is equivalent to increasing the field at fixed dis- 
order. At low field or low disorder, the ground state is 
a Bragg glass with 9 ~ -l-l.O, exactly as with screening. 
Without screening, the main difference is that the high 
field, large disorder, phase is a true vortex glass with 
stiffness exponent 9 ~ -I-0.30, as proposed by Fisher et 
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We tentatively conclude that a true superconduct- 
ing vortex glass phase does not exist in except in the 
absence of screening (A oo). Our understanding of the 
experimental consequences for real systems with meso- 
scopic penetration depths A ~ C'(10'^)A is lacking and 



may be of some interest. 

Computations were performed at the Theoretical 
Physics Computing Facility at Brown University. JMK 
thanks A. Vallat and B. Grossman for many discussions 
about spin glasses, best twists etc. 
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